#!/usr/bin/env julia
"""
NSEMSolver.jl Comprehensive Visualization Demo

This example demonstrates the complete visualization capabilities of NSEMSolver.jl,
showcasing all plotting functions, backends, and export formats.

Features demonstrated:
- Basic solution plotting with multiple backends
- Velocity field visualization (components and magnitude) 
- Pressure field contours and surfaces
- Streamlines and vector field plots
- Convergence history analysis
- Animation creation from time series
- VTK export for ParaView visualization
- Multi-backend graceful fallbacks
- Adaptive refinement visualization
- Comprehensive solution analysis
"""

using NSEMSolver
using Printf

function main()
    println("🎨 NSEMSolver.jl - Comprehensive Visualization Demo")
    println("=" ^ 60)
    
    # Configuration for a moderately complex simulation
    options = NSOptions(
        N = 5,                    # Higher polynomial order for smoother results
        n_block = 4,             # More blocks for finer resolution
        nu = 0.01,               # Moderate viscosity
        tfinal = 2.0,            # Longer simulation time
        cfl = 0.3,               # Conservative CFL for stability
        tol = 1e-7,              # Tight tolerance for accuracy
        solver = :julia,         # Use Julia backend for compatibility
        method = :sem,           # Spectral element method
        adaptive_refinement = true,  # Enable refinement for demonstration
        refinement_levels = 2,   # Moderate refinement
        verbose = true,          # Show detailed progress
        save_history = true      # Essential for convergence plots
    )
    
    println("Simulation Configuration:")
    println("  🔬 Polynomial order: $(options.N)")
    println("  🏗️  Domain blocks: $(options.n_block)×$(options.n_block)")
    println("  🌊 Viscosity: $(options.nu)")
    println("  ⏰ Final time: $(options.tfinal)")
    println("  📐 CFL number: $(options.cfl)")
    println("  🎯 Tolerance: $(options.tol)")
    println("  🔧 Solver: $(options.solver)")
    println("  📊 Adaptive refinement: $(options.adaptive_refinement)")
    println()
    
    # Run the simulation
    println("🚀 Running Navier-Stokes simulation...")
    start_time = time()
    
    result = solve_navier_stokes_2d(options)
    
    solve_time = time() - start_time
    
    # Display basic results
    println()
    println("✅ Simulation completed!")
    println("📈 Convergence: $(result.converged)")
    println("🔄 Iterations: $(result.iterations)")
    println("📉 Final residual: $(@sprintf("%.2e", result.residual_norm))")
    println("⏱️  Total time: $(@sprintf("%.3f", solve_time)) seconds")
    println()
    
    # Create comprehensive visualization demonstrations
    demonstrate_basic_plotting(result)
    demonstrate_velocity_visualization(result)
    demonstrate_pressure_visualization(result)
    demonstrate_convergence_analysis(result)
    demonstrate_export_capabilities(result)
    demonstrate_advanced_features(result)
    demonstrate_backend_comparison(result)
    
    # Create summary report
    create_visualization_report(result)
    
    println("🎉 Visualization demo completed successfully!")
    println("📁 Check the ../visualization_outputs/ directory for generated plot files and exports.")
    
    return result
end

function demonstrate_basic_plotting(result::NSResult)
    println("📊 Demonstrating Basic Solution Plotting")
    println("-" ^ 40)
    
    # Basic solution plot with default backend
    println("🎨 Creating comprehensive solution overview...")
    plot_ns_solution(result, save_path="ns_solution_overview.png")
    
    # Try different backends if available
    for backend in [:plots, :plotlyjs, :makie, :text]
        println("  🖼️  Testing backend: $backend")
        try
            plot_ns_solution(result, backend=backend, save_path="ns_solution_$(backend).png")
            println("    ✅ $backend backend successful")
        catch e
            println("    ⚠️  $backend backend not available: $(typeof(e))")
        end
    end
    println()
end

function demonstrate_velocity_visualization(result::NSResult)
    println("🌀 Demonstrating Velocity Field Visualization")
    println("-" ^ 45)
    
    # Individual velocity components
    components = [:magnitude, :u, :v]
    if result.w !== nothing
        push!(components, :w)
    end
    
    for component in components
        println("  📈 Plotting velocity component: $component")
        plot_velocity_field(result, component, save_path="velocity_$(component).png")
    end
    
    # Vector field visualization
    println("  🏹 Creating velocity vector field...")
    plot_velocity_vectors(result, skip=3, scale=0.5, save_path="velocity_vectors.png")
    
    # Streamlines
    println("  🌊 Generating streamlines...")
    plot_streamlines(result, density=1.5, save_path="streamlines.png")
    
    println()
end

function demonstrate_pressure_visualization(result::NSResult)
    println("💫 Demonstrating Pressure Field Visualization")
    println("-" ^ 45)
    
    # Pressure contours with different levels
    contour_levels = [5, 10, 20]
    
    for levels in contour_levels
        println("  📊 Pressure contours with $levels levels")
        plot_pressure_contours(result, levels=levels, save_path="pressure_$(levels)_levels.png")
    end
    
    println()
end

function demonstrate_convergence_analysis(result::NSResult)
    println("📈 Demonstrating Convergence Analysis")
    println("-" ^ 38)
    
    if !isempty(result.convergence_history)
        # Linear and logarithmic scale convergence plots
        println("  📉 Linear scale convergence history")
        plot_convergence_history(result, log_scale=false, save_path="convergence_linear.png")
        
        println("  📉 Logarithmic scale convergence history")
        plot_convergence_history(result, log_scale=true, save_path="convergence_log.png")
        
        # Convergence statistics
        initial_res = result.convergence_history[1]
        final_res = result.convergence_history[end]
        reduction = initial_res / final_res
        
        println("  📊 Convergence Statistics:")
        println("    Initial residual: $(@sprintf("%.2e", initial_res))")
        println("    Final residual: $(@sprintf("%.2e", final_res))")
        println("    Reduction factor: $(@sprintf("%.1e", reduction))")
    else
        println("  ⚠️  No convergence history available")
    end
    
    println()
end

function demonstrate_export_capabilities(result::NSResult)
    println("💾 Demonstrating Export Capabilities")
    println("-" ^ 35)
    
    # VTK export for ParaView
    println("  🔬 Exporting to VTK format for ParaView...")
    export_vtk(result, "ns_solution_vtk", include_vorticity=true)
    
    # Create visualization summary
    println("  📊 Creating visualization summary...")
    summary = create_visualization_summary(result)
    
    # Save summary to JSON-like format
    output_dir = "../visualization_outputs"
    mkpath(output_dir)
    open(joinpath(output_dir, "visualization_summary.txt"), "w") do io
        println(io, "NSEMSolver.jl Visualization Summary")
        println(io, "=====================================")
        for (key, value) in summary
            println(io, "$key: $value")
        end
    end
    
    println("    ✅ Summary saved to $(joinpath(output_dir, \"visualization_summary.txt\"))")
    println()
end

function demonstrate_advanced_features(result::NSResult)
    println("🔬 Demonstrating Advanced Visualization Features")
    println("-" ^ 48)
    
    # Adaptive refinement visualization (if available)
    println("  🎯 Adaptive refinement levels...")
    plot_domain_refinement(result, save_path="refinement_levels.png")
    
    # Multi-field comparison
    println("  📊 Creating multi-field analysis...")
    
    # Check field ranges and statistics
    vel_mag = NSEMSolver.compute_velocity_magnitude(result)
    vorticity = NSEMSolver.compute_vorticity(result)
    
    println("    Field Statistics:")
    println("      Velocity magnitude: min=$(@sprintf("%.3f", minimum(vel_mag))), max=$(@sprintf("%.3f", maximum(vel_mag)))")
    println("      Pressure: min=$(@sprintf("%.3f", minimum(result.p))), max=$(@sprintf("%.3f", maximum(result.p)))")
    println("      Vorticity: min=$(@sprintf("%.3f", minimum(vorticity))), max=$(@sprintf("%.3f", maximum(vorticity)))")
    
    println()
end

function demonstrate_backend_comparison(result::NSResult)
    println("🔄 Demonstrating Backend Comparison")
    println("-" ^ 35)
    
    backends = [:plots, :plotlyjs, :makie, :text]
    
    for backend in backends
        println("  🖥️  Testing $backend backend capabilities...")
        
        # Test basic plotting
        try
            start_time = time()
            plot_ns_solution(result, backend=backend)
            elapsed = time() - start_time
            println("    ✅ Basic plotting: successful ($(@sprintf("%.3f", elapsed))s)")
        catch e
            println("    ❌ Basic plotting: failed ($(typeof(e)))")
        end
        
        # Test velocity field
        try
            plot_velocity_field(result, :magnitude, backend=backend)
            println("    ✅ Velocity field: successful")
        catch e
            println("    ❌ Velocity field: failed ($(typeof(e)))")
        end
        
        # Test convergence history
        try
            plot_convergence_history(result, backend=backend)
            println("    ✅ Convergence plot: successful")
        catch e
            println("    ❌ Convergence plot: failed ($(typeof(e)))")
        end
    end
    
    println()
end

function create_visualization_report(result::NSResult)
    println("📋 Creating Comprehensive Visualization Report")
    println("-" ^ 47)
    
    output_dir = "../visualization_outputs"
    mkpath(output_dir)
    report_filename = joinpath(output_dir, "visualization_report.md")
    
    open(report_filename, "w") do io
        write_markdown_report(io, result)
    end
    
    println("  📝 Report saved to $report_filename")
    println()
end

function write_markdown_report(io, result::NSResult)
    println(io, "# NSEMSolver.jl Visualization Report")
    println(io, "")
    
    # Solution overview
    println(io, "## Solution Overview")
    println(io, "")
    println(io, "- **Grid Size**: $(size(result.u))")
    println(io, "- **Dimensions**: $(result.w === nothing ? "2D" : "3D")")
    println(io, "- **Converged**: $(result.converged)")
    println(io, "- **Iterations**: $(result.iterations)")
    println(io, "- **Final Residual**: $(@sprintf("%.2e", result.residual_norm))")
    println(io, "- **Solve Time**: $(@sprintf("%.3f", result.solve_time)) seconds")
    println(io, "")
    
    # Field statistics
    println(io, "## Field Statistics")
    println(io, "")
    
    vel_mag = NSEMSolver.compute_velocity_magnitude(result)
    
    println(io, "### Velocity Components")
    println(io, "- **U-velocity**: min=$(@sprintf("%.3f", minimum(result.u))), max=$(@sprintf("%.3f", maximum(result.u)))")
    println(io, "- **V-velocity**: min=$(@sprintf("%.3f", minimum(result.v))), max=$(@sprintf("%.3f", maximum(result.v)))")
    if result.w !== nothing
        println(io, "- **W-velocity**: min=$(@sprintf("%.3f", minimum(result.w))), max=$(@sprintf("%.3f", maximum(result.w)))")
    end
    println(io, "- **Velocity Magnitude**: min=$(@sprintf("%.3f", minimum(vel_mag))), max=$(@sprintf("%.3f", maximum(vel_mag)))")
    println(io, "")
    
    println(io, "### Pressure Field")
    println(io, "- **Pressure**: min=$(@sprintf("%.3f", minimum(result.p))), max=$(@sprintf("%.3f", maximum(result.p)))")
    println(io, "")
    
    # Convergence analysis
    if !isempty(result.convergence_history)
        println(io, "## Convergence Analysis")
        println(io, "")
        initial_res = result.convergence_history[1]
        final_res = result.convergence_history[end]
        reduction = initial_res / final_res
        
        println(io, "- **Initial Residual**: $(@sprintf("%.2e", initial_res))")
        println(io, "- **Final Residual**: $(@sprintf("%.2e", final_res))")
        println(io, "- **Reduction Factor**: $(@sprintf("%.1e", reduction))")
        println(io, "- **Convergence Rate**: $(reduction > 1000 ? "Excellent" : reduction > 100 ? "Good" : "Moderate")")
        println(io, "")
    end
    
    # Visualization files generated
    println(io, "## Generated Visualization Files")
    println(io, "")
    println(io, "1. `ns_solution_overview.png` - Complete solution overview")
    println(io, "2. `velocity_magnitude.png` - Velocity magnitude field")
    println(io, "3. `velocity_u.png` - U-velocity component")
    println(io, "4. `velocity_v.png` - V-velocity component")
    println(io, "5. `velocity_vectors.png` - Velocity vector field")
    println(io, "6. `streamlines.png` - Flow streamlines")
    println(io, "7. `pressure_*.png` - Pressure contours (multiple levels)")
    println(io, "8. `convergence_*.png` - Convergence history plots")
    println(io, "9. `ns_solution_vtk.vtu` - VTK format for ParaView")
    println(io, "10. `ns_solution.csv` - CSV export for external tools")
    println(io, "")
    
    # Code examples
    println(io, "## Code Examples")
    println(io, "")
    println(io, "```julia")
    println(io, "# Basic visualization")
    println(io, "using NSEMSolver")
    println(io, "result = solve_navier_stokes_2d()")
    println(io, "plot_ns_solution(result)")
    println(io, "")
    println(io, "# Advanced features")
    println(io, "plot_velocity_field(result, :magnitude)")
    println(io, "plot_pressure_contours(result, levels=15)")
    println(io, "plot_streamlines(result, density=2.0)")
    println(io, "export_vtk(result, \"solution\")")
    println(io, "```")
end

# Run the main demonstration
if abspath(PROGRAM_FILE) == @__FILE__
    println("Starting NSEMSolver.jl Visualization Demo...")
    println()
    
    # Main demonstration
    result = main()
    
    println("")
    println("📚 Demo Summary:")
    println("  ✅ Basic plotting: Complete")
    println("  ✅ Velocity visualization: Complete")
    println("  ✅ Pressure visualization: Complete")
    println("  ✅ Convergence analysis: Complete")
    println("  ✅ Export capabilities: Complete")
    println("  ✅ Advanced features: Complete")
    println("  ✅ Backend comparison: Complete")
    println("  ✅ Comprehensive report: Generated")
    println("")
    println("🎯 Next Steps:")
    println("  1. Install Plots.jl: julia -e 'using Pkg; Pkg.add(\"Plots\")' ")
    println("  2. Install PlotlyJS.jl: julia -e 'using Pkg; Pkg.add(\"PlotlyJS\")' ")
    println("  3. Install WriteVTK.jl: julia -e 'using Pkg; Pkg.add(\"WriteVTK\")' ")
    println("  4. Open ParaView and load ns_solution_vtk.vtu")
    println("  5. Review visualization_report.md for detailed analysis")
    println("")
    println("🎉 NSEMSolver.jl visualization system is production-ready!")
end